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On  Properties  of  Contours  of 
lYilinear  Scalar  Fields 


Holger  Theisel 


Abstract.  We  study  properties  of  contour  surfaces  of  trilinear  scalar 
fields,  and  give  a classification  based  on  how  many  unconnected  surface 
parts  they  consist  of.  Furthermore,  we  introduce  the  concept  of  the  seg- 
ment number  of  a voxel.  The  segment  number  is  a threshold-independent 
measure  which  estimates  how  complicated  the  contours  inside  the  voxel 
are  expected  to  be.  Finally,  we  give  necessary  and  sufficient  conditions  for 
a voxel  to  have  a segment  number  of  1.  These  conditions  are  applied  to 
analyze  a computer  tomography  data  set. 


§1.  Introduction 

Contours  (isosurfaces)  of  trilinear  scalar  fields  are  treated  in  a variety  of  appli- 
cations. For  instance,  the  data  used  in  volume  visualization  usually  consists 
of  a number  of  scalars  defined  at  certain  grid  points;  between  the  grid  points 
a piecewise  trilinear  interpolation  of  the  scalar  field  is  applied. 

Given  a voxel  V = [0,  l]3,  the  trilinear  scalar  field  is  defined  by  setting 
the  values  Cijk(i,j,  k € {0, 1})  of  the  field  at  the  corners  of  V.  Then  the  scalar 
field  is  defined  as 

s(u,  V,  w)  — (1  - u)  ■ (1  - v)  ■ (1  - w)  ■ Cooo  + (1  - u)  ■ (1  - v)  ■ w ■ c0oi 

+ (1  - It)  • V ■ (1  - w)  ■ C010  + (1  — u)  ■ V ■ W ■ Coil  ^ 

+ u ■ (1  - v)  ■ (1  - w)  ■ ci oo  + u ■ (1  - v)  ■ w ■ c10i 
+ u • v • (1  — w)  ■ cuo  + u ■ v ■ w ■ cm- 


Figure  la  illustrates  this.  A contour  of  V is  defined  by  s(u,v,w)  = r =const 
for  a certain  threshold  r.  Figure  lb  shows  an  example  of  a contour  of  (1). 

There  are  a number  of  algorithms  to  produce  a triangular  approximation 
of  a contour  of  (1).  Of  these,  the  Marching  Cubes  (MC)  method  ([3]  and  [4]) 
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Fig.  1.  a)  Voxel  V;  b)  a contour  in  V;  c)  result  of  MC. 

is  the  most  popular.  Figure  lc  shows  the  resulting  triangular  approximation 
of  the  contour  shown  in  Figure  lb  using  the  Marching  Cubes  method. 

The  Marching  Cubes  algorithm  distinguishes  several  cases  where  some  of 
them  are  harder  to  treat  than  others.  In  this  paper  we  introduce  a measure 
of  how  costly  in  terms  of  computing  time  the  MC  algorithm  inside  a certain 
voxel  is  expected  to  be.  This  characterization  of  a voxel  - called  segment 
number  - is  independent  of  a particular  threshold.  It  estimates  the  costs  of  a 
Marching  Cubes  algorithm  for  varying  thresholds. 

As  already  stated  in  [2],  the  contour  of  (1)  is  a rational  cubic  surface.  In 
[2]  this  surface  is  approximated  by  a collection  of  rational  quadratic  triangular 
patches. 

Section  2 of  this  paper  studies  the  contours  of  (1)  in  the  domain  1R3. 
We  give  a classification  based  on  how  many  unconnected  surface  parts  the 
contours  consist  of.  Sections  3 and  4 focus  on  contours  of  (1)  inside  a certain 
voxel.  Section  3 introduces  the  concept  of  segment  number  as  a measure  of 
how  simply  a voxel  can  be  treated  by  an  MC  algorithm.  In  Section  4,  necessary 
and  sufficient  geometric  conditions  for  a voxel  to  have  a segment  number  of 
1 are  shown.  In  Section  5,  the  number  of  voxels  with  a segment  number  of  1 
are  computed  for  a real  volume  data  set. 

§2.  Classification  of  the  Contour  in  R3 

In  this  section  we  consider  the  contour  of  (1)  not  in  a particular  voxel  but  in 
the  domain  IR3 . In  general,  the  contour  consists  of  a number  of  surface  parts 
which  are  not  connected  to  each  other.  Before  we  classify  the  contours  of  (1) 
by  the  number  of  unconnected  surface  parts,  we  apply  a translation  of  the 
coordinate  system  as  shown  in  Figure  2.  Choosing 


P = cooi  + coio  + cioo  + cm  — cooo 
( Cooo  + Coil  — Cool  — Coio  \ 

1 


Coil  — Cioi  — Clio 


Po  = -- 
p 


Cooo  + Cioi  — Cool  — Cioo 
\cq00  + Clio  — Coio  — Cioo  ) 


we  obtain  for  (1) 


s = a-  u + b-  v + c-  w + d-  u-  v-  w + e 


(2) 
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Fig.  2.  Translating  the  coordinate  system  of  a voxel. 

with 

(Cm  - C0ii)  • (cioo  - Cooo)  - (cno  - c0io)  • (cioi  - Cool) 

a = 

V 

(cm  - Cioi)  • (coio  - Cooo)  - (cno  - Cioo)  • (coil  - Cool) 

o = — 

p 

(cill  - Clio)  • (cool  - Cooo)  - (cioi  - Cioo)  • (coil  - Coio) 

c = 

p 

d — Pi 

where  e is  a certain  constant.  Thus,  we  only  have  to  analyze 

s(u,  v,w)  = a-  u + b-  v + c-  w + d-  u-  v-  w = r = const  (3) 

in  1R3.  A classification  of  (3)  can  be  achieved  by  rewriting  (3)  as  w = r~Xd-ZsT 
and  comparing  the  zeros  of  the  numerator  and  denominator  function.  The 
zeros  of  the  numerator  function  form  a line  in  the  u — v— plane,  whereas  the 
zeros  of  the  denominator  function  give  a hyperbola.  Studying  their  interplay 
gives  the  following  classification: 

case  1:  abed  < 0,  d ^ 0 : 

case  1.1:  r2  > -^p:  (3)  gives  3 unconnected  surface  parts 
case  1.2:  r2  < -pp:  (3)  gives  2 unconnected  surface  parts 
case  2:  abed  < 0:  (3)  consists  of  1 connected  part 
case  3:  abed  = 0,  d ^ 0: 
case  3.1:  r^0: 

case  3.1.1:  06  / 0,c  = 0:  (3)  gives  2 unconnected  surface  parts 
case  3.1.2:  a ^ 0, 6 = c = 0:  (3)  gives  3 unconnected  surface  parts 

case  3.1.3:  a = b = c = 0:  (3)  gives  4 unconnected  surface  parts 

case  3.2:  r = 0 : 

case  3.2.1:  a6  ^ 0,c  = 0:  (3)  gives  3 unconnected  surface  parts 
case  3.2.2:  a / 0, 6 = c = 0:  (3)  gives  3 parts  intersecting  each  other 

case  3.2.3:  a = b = c = 0:  (3)  gives  3 perpendicular  planes. 


Figure  3 illustrates  these  cases. 
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Fig.  3.  Classification  of  the  contours  of  (3)  in  fi'!  . 

§3.  Segment  Number  of  a Voxel 

We  now  study  the  contour  of  (3)  in  a particular  voxel  V = [tig,Uo  + 1]  x 
[uo,vo  + 1]  x [wo,wo  + 1].  Unfortunately,  the  results  of  Section  2 are  not 
directly  applicable  here  because  one  connected  surface  part  may  intersect  V 
more  than  once. 

Varying  the  threshold  r in  (3),  the  contours  change.  So  does  the  number 
of  unconnected  surface  parts  of  the  contour. 

Definition  1.  Given  the  triline  ar  scalar  field  s(u,  v,w)  = a-u  + b-  v + c-  w + 
d-u-v-w  in  the  domain  of  the  voxel  V = [uo,  u0  + 1]  x [v0,  fo  + 1]  x [wo,  wo  + 1], 
the  segment  number  S(V ) of  V is  the  maximal  number  of  unconnected  surface 
parts  of  the  contour  s(u,  v,w)  = r =const  in  V for  any  threshold  r. 

Figure  4 gives  an  example  of  a voxel  V with  S(V)  = 1.  Increasing  the 
value  of  r,  the  isosurface  ’’moves”  through  the  voxel.  It  consists  of  at  most 
one  connected  part  for  any  r.  Figure  5 shows  a voxel  with  S(V)  = 4.  Here 
the  contours  consist  of  up  to  4 unconnected  parts. 

The  segment  number  is  a threshold-independent  characterization  of  a 
voxel  V.  For  any  V we  get  5(V)  € {1,2,3, 4}.  For  visualization  purposes, 
voxels  with  a segment  number  1 are  of  particular  interest.  As  shown  in  the 
example  of  Figure  4,  they  have  a nice  behavior  while  varying  r.  In  fact,  for 


On  Contours  of  Trilinear  Scalar  Fields 


407 


Fig.  5.  Contours  of  a voxel  with  S(V)  = 4. 


any  r the  contour  consists  of  only  one  connected  surface  part  inside  V.  Thus, 
accelerated  Marching  Cubes  methods  may  apply  to  them.  Moreover,  adjacent 
voxels  with  S(V)  = 1 may  be  merged  to  form  one  bigger  voxel  before  applying 
Marching  Cubes  methods.  So  it  makes  sense  to  search  for  geometric  conditions 
for  a voxel  V to  have  S(V)  = 1. 


§4.  Geometric  Conditions  for  S(V)  = 1 

In  this  section  we  give  necessary  and  sufficient  geometric  conditions  for  a 
voxel  to  have  S(V)  = 1.  Again,  we  consider  the  contour  of  (3)  in  the  voxel 
V = [u0,Mo  + 1]  X [vo,UO  + 1]  x [w0,  w0  + 1]. 

To  formulate  the  conditions  for  S(V)  = 1,  we  need  to  introduce  the 
concept  of  characteristic  hyperbolas.  The  first  characteristic  hyperbola  hi  in 
]R3  is  defined  by  the  condition  svw(u,v,w ) = 0 in  (3).  h\  can  be  written  as  a 
rational  quadratic  Bezier  curve  described  by  two  control  vectors  6q,  b\  and  a 
control  point  b\  (see  [1]).  For  hi  we  obtain 


f{-Abc)/d\ 

( °\ 

bl  = 

0 

, *1  = 

0 , b\  = 

1/6 

l o ) 

w 

w«/ 
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Fig.  6.  Location  of  characteristic  hyperbolas;  a),b):  abed  < 0;  c),d):  abed  > 0. 
where  w\  is  the  weight  of  6} . Then  we  obtain 


hi(t) 


blBl{t)  + w\b\Bl{t)  + b\Bl{t) 
w\Bl(t) 


In  a similar  way  we  define  the  characteristic  hyperbola  h2  by  suw(u,  v,  w)  = 0, 
and  h3  by  su„(«,  v,  w ) = 0.  The  Bezier  point  6j  with  the  corresponding  weight 
w\  and  the  control  vectors  6„ , b?2  describing  h2  are 


h3  is  described  by 


( 


\{-4ab)/d/ 


f°\ 

( l!a\ 

*1=0,  bl  = 

0 

w 

\l/‘) 

II 

CO  -H 

o o 

l3 

5 ^2  — 

^ l/a\ 
1/6 

l 0 / 

vi\  = l. 


w?  = 1. 


If  a • b ■ c ■ d < 0 then  hi,  h2,  h3  intersect  in  two  common  points.  Figures  6a 
and  6b  illustrate  this  situation  from  two  different  viewpoints.  If  a ■ b ■ c • d > 0, 
then  hi,h2,h3  do  not  have  any  intersections.  Figures  6c  and  6d  show  this 
from  different  viewpoints.  The  degenerate  case  a ■ b ■ c ■ d — 0 is  omitted  here. 

To  formulate  conditions  for  S(V)  = 1,  we  have  to  classify  the  faces  of  V. 
Given  the  voxel  V = [«0,  u0  + l]  x [v0,  «o  + l]  x [wo,  u>0  + l],  let  — {( u,v,w ) G 

V : u = Mo  V u = uq  + 1},  f2  = {(it,v,w)  G V : v = v0  V v = vq  + 1},  and 
/ 3 = {(u,v,w)  G V : w = wq  V w = wg  + 1}.  See  Figure  7 for  an  illustration 
of  the  faces. 


Theorem  1.  Let  V = [moj^o  + 1]  x [vo,uo  + 1]  x [wo,wo  + 1]  be  a voxel  in 
the  scalar  held  defined  by  (3).  Then  the  condition  S(V)  = 1 is  equivalent  to 
the  three  conditions  hi  PI  f x = 0 and  h2  fl  f2  = 0 and  h3  fl  f3  = 0. 

Figure  8 illustrates  the  idea  of  the  proof.  Suppose  h3  intersects  f 3 as 
shown  in  Figure  8a.  Figure  8b  is  a magnification  of  the  voxel  and  h3  in 
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Fig.  8.  Proof  idea  of  Theorem  1. 

Figure  8a.  We  compute  the  intersection  point  of  A3  and  /3,  and  consider 
the  contour  passing  through  this  point.  As  shown  in  Figure  8a,  this  contour 
consists  of  at  least  two  surface  parts. 

For  the  proof  of  the  converse  statement  of  Theorem  1,  we  assume  that  for 
a certain  threshold  r the  contour  consists  of  at  least  two  unconnected  surface 
parts.  Then  we  can  find  a face  of  F which  has  two  intersection  curves  with  the 
contour.  (In  the  worst  case  we  have  to  vary  r to  find  such  a face).  (Figure  8c 
shows  two  surface  parts  of  the  contour  which  produce  two  intersection  curves 
in  the  upper  face  of  /3).  Then  we  can  find  a point  on  this  face  which  is  the 
intersection  point  with  the  corresponding  characteristic  hyperbola.  (In  Figure 
8c,  the  marked  point  on  the  upper  part  of  /3  is  the  intersection  with  A3). 


§5.  Results  and  Future  Work 

We  have  tested  the  voxels  of  a CT  test  data  set  for  the  property  5(F)  = 1. 
The  data  set  consists  of  255  x 255  x 108  = 7, 022, 700  voxels.  Figure  9 shows 
a slice  through  the  data  set. 

In  the  raw  data  we  found  1,978,711  voxels  with  5(F)  = 1 (28  %).  After 
some  noise  reducing  filter  operations  on  the  data,  we  detected  4,833,063  voxels 
with  5(F)  = 1 (69  %).  This  shows  that  there  is  a reasonable  number  of  voxels 
with  5(F)  = 1 to  pay  special  attention  to  them. 

In  the  future  we  plan  to  develop  algorithms  to  merge  voxels  with  5(F)  = 1 
to  form  bigger  voxels  before  starting  the  Marching  Cubes  algorithm. 
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Fig.  9.  Slice  through  the  test  data  set. 
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